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I. INTRODUCTION 

The electron gas with long-range Coulomb repulsion 
and the Hubbard model with screened on-site repulsion 
are two widely studied models that each describe large 
classes of materials. Both models are generally studied 
with quite different theoretical methods. However, by 
gradually increasing the range of the interaction in the 
Hubbard model and by reducing the density to average 
out lattice effects, one should be able to go continuously 
from the Hubbard to the Coulomb gas model. Is there 
one unified theoretical framework that allows us to treat 
both limiting cases? 

To begin to answer this question, we generalize 
the Two-Particle Self-Consistent approach (TPSC)±^ 
to study the extended Hubbard model, that includes 
nearest-neighbor repulsion in addition to the usual on- 
site repulsion. This model is interesting in its own 
right, independently of the above-mentioned theoreti- 
cal question. Indeed, the extended Hubbard model al- 
lows one to study materials where competition between 
charge and spin order manifest themselves. In high- 
temperature superconductors, where screening is not per- 
fect, understanding the extended Hubbard model is also 
of paramount importance. 

Let us first motivate further our focus on the TPSC 
approach and then come back to the interesting physi- 
cal phenomena that manifest themselves in the extended 
Hubbard model. Judging from comparisons with bench- 
mark Quantum Monte Carlo calculationsji*2i2i4i£i& (in the 
absence of exact solutions), the TPSC approach pro- 
vides us with the most accurate approximate solution 
to the Hubbard model at weak to intermediate coupling. 
A detailed critical comparison with other methods such 
as the Random Phase Approximation (RPA), the self- 
consistent renormalized theory and the fluctuation ex- 
change approximation is given in Ref. |2). In particular, 
TPSC satisfies the Pauli principle, the Mcrmin- Wagner 
theorem in two-dimensions includes Kanamori-BriickneH' 
screening and is non-perturbative although limited to in- 
teraction strengths less than the bandwidth. The TPSC 



approach is a close relative to the Singwi, Tosi, Land, 
Sjolander (STLS) method^ used in the electron gas prob- 
lem. The STLS approximation has been first introduced 
to describe the structure functions of an electron liquid 
where it provided much better results compared to the 
RPA. This approximation has been applied to a vari- 
ety of systems that contain Fermions, Bosons or classical 
particles in all physical dimensions and different geome- 
tries. Starting with the equation of motion for one-body 
density operator, the authors were faced with the well- 
known problem that the two-body density operator ap- 
peared in their equation. They solved the problem ap- 
proximately by replacing the two-body density operator 
by a product of two one-body density operators and then 
correcting the result with the pair correlation function. 
The result of this factorization appears as a correction 
in the response functions of the system, the so-called lo- 
cal field factor. This factor is then determined by using 
a sum-rule derived from the fluctuation-dissipation the- 
orem. The present paper will give a new point of view 
on the STLS method by comparing it to the TPSC ap- 
proach. As we will show in more details, the main dif- 
ference between the latter method and the STLS one, is 
the way we factorize the two-body density operator. It 
seems that for local vs non-local potentials, it is more 
accurate to use, respectively, local or non-local factoriza- 
tion and, as we will show in detail, a local factorization 
leads to better results for the extended Hubbard model. 
That is not all. It will become clear in the formalism 
used to derive TPSC that the STLS method also ne- 
glects some higher-order correlation functions. The same 
type of approximation will be necessary to be able to 
generalize TPSC to treat the extended Hubbard models. 
Otherwise, as in STLS, there is a shortage of sum-rules 
or conditions to find all unknowns that appear in the 
method. 

We believe that accuracy of the approximation is cru- 
cial for a real understanding of physical properties and 
for meaningful comparisons to experiments. Bad ap- 
proximations that agree with experiment only lead us 
astray. That is why we will benchmark our extension of 
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the TPSC approach against the highly accurate results 
that can obtained by Quantum Monte Carlo (QMC) sim- 
ulations. We will carefully analyze the approximations 
involved in the method and discuss other possibilities for 
improvement. 

Back to the extended Hubbard model. This model 
has a long history, so we can only discuss a small sam- 
ple of the relevant literature. At half-filling, when the 
on-site interaction strength U tends to infinity so that 
super-exchange At 2 /U vanishes, the effect of the nearest- 
neighbor repulsion V is to lead to an effective ferro- 
magnetic interaction between localized spins. This is 
the physics of the so-called direct-exchange mechanism. 
The physics is quite different when U and AV are of 
the same order and both in the weak to intermediate- 
coupling regime, namely less than or of the order of the 
bandwidth (W = 8t in d = 2). In that case, there 
is a competition between staggered charge and spin or- 
ders. That charge ordering phenomenon is particularly 
relevant for manganites, vanadates and various organic 
conductors, as discussed in a recent theoretical paper- 
that uses a new correlator-projection method. The rele- 
vant theoretical literature for these compounds focuses on 
the square latticeAi*i2iiii and on ladders for the quarter- 
filled cas o 14 ! 15 . The competition between charge and spin 
orders has also been studied in one-^iLi^ two-^2i22*2I, 
three— and higher dimensions^ at various fillings. The 
combined effect of charge fluctuations in addition to the 
usual spin fluctuation in favouring one type or another of 
unconventional superconductivity has also been studied 



using this model i 25 i 26 . 

In the present paper, we are interested in the possibil- 
ity that a generalization of TPSC to the extended Hub- 
bard model can lead us to accurate estimations of the 
charge and spin structure factors and susceptibilities at 
finite temperature. We will use the QMC calculations 
of Ref. (f30r) as a benchmark. We note that methods 
that have been quite successful in one- or infinite dimen- 
sion are generally not applicable in the two-dimensional 
case that we will consider. In d = 2, continuous sym- 
metries can be broken only at zero temperature and, in 
addition, wave-vector dependencies that are neglected in 
high-dimension are generally not negligible. 

In the following, we first present the theory and give 
the details of the calculation based on the functional 
derivative of the Dyson equation which gives us the re- 
sponse functions of the system. We also provide the equa- 
tion of motion for the Wigner distribution function to 
show that the two different methods basically lead to the 
same set of equations. This also allows us to discuss the 
different types of factorizations. In Sec. lIIII we present the 
results of our numerical calculations and compare them 
with QMC results to find the region where the method 
works properly or is precise enough. Finally, we discuss 
the influence of nearest-neighbor interaction V on spin 
and charge fluctuations. 

II. THEORY 

We first introduce the extended Hubbard Hamiltonian, 



H = 



) + U Yl ni T n H + v E n ™ n . 



E 



"J<r 



(1) 



where c- lcr (el ) are annihilation (creation) operator for 
electrons of spin a at site i, n\ a is the density operator, 
and t is the hopping matrix element. The quantities U 
and V are the on-site and nearest-neighbor interactions 
respectively and fj, is the chemical potential. Although we 
restrict ourselves to nearest-neighbor hopping, the gen- 
eralization to an arbitrary hopping matrix tjj will be ob- 
vious to the reader. It only modifies the non-interacting 
dispersion relation. One can generalize the formalism to 
a system with longer interaction terms and also to a sys- 
tem with many bands. In the following, we first derive 
the TPSC approach using functional derivatives^ and 
then return to the approach that is more usual with the 
STLS approximation 8 , namely the equation of motion 
for the one-body Wigner distribution function. These 
two derivations allow a deeper insight into the nature 
of the approximations. The reader may also choose the 
approach he is more familiar with. 



A. Functional derivative approach 

Following functional methods of the Schwinger 
schooi22i2iiM, we begin with the generating function 
lnZ^o-] with source fields <p a in the grand canonical 
ensemble 

Z\$ a \ = -Tr [e-f )H T T S{p)\ (2) 

where f3 = 1/T, T is the temperature while S is defined 
as follows 

]nS(/3) = -J2 dTdr' C L(T) Cja (T')^(iJ,r,r') (3) 

y,> 

where the brackets () represents a thermal average in 
the canonical ensemble, T T is the time-ordering operator, 
and t is the imaginary time. The Green function can 
be calculated from the first functional derivative of the 
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generating function Z [(/>„] as follows 

SlnZfa] _ Tr [ e -^T T 5(/?)4(2) C(7 (l)] 



G CT (1,2) 



M2,l) 



Tr [e-P B T T S(0)] 



(4) 

where we have introduced the short-hand 1 to stand for 
both the site position and the corresponding imaginary 



J 



time, as in the equation 

G CT (1,2) = - (TrClM&rfo)) = - (T T c CT (l)4(2)) . 

(5) 

The equation of motion for the Green function has the 
following forr" 2 ' 31 ' 33 



G^(l, 2)G CT (2, 3) = 5(1, 3) + £ CT (1, 2)G CT (2, 3) + ^(1, 2)G CT (2, 3) (6) 

where G " 1 (l,2) = [J^ — e(l, 2)] J(l , 2) is the non-interacting Green function, S CT is the self-energy and the bar is a 
short-hand for J dr. The above equation is nothing more than the Dyson equation that can be obtained by the 
diagrammatic technique. It can also be written in the form 



G- 1 (1,2)=G 1 (1,2)-^(1,2)-S CT (1,2), 



(7) 



with the self-energy 

£ CT (1,2) = -[/(T T 4(l)c*(l) C(T (l)4(3)) G-\l,2) - VJ2 (^4,(1 + o)<v(l + a)e CT (l)4(3)) 0^(3,2) (8) 



where a = — a and the summation on a runs over the nearest-neighbor sites of the site 1 (the imaginary time is the 
same at 1 + a as at 1). One needs an approximation to deal with the above two-body density operator. By analogy 
with the factorization pioneered by Singwi et. al £ we write 

S CT (1, 2) = UG a (l, 1 + )G CT (1, 3)G; 1 (3, 2) 5 ^(1, 1) + G„, (1 + a, 1 + a+)G (T (l, 3)G; X (3, 2)<w(l, 1 + «) 

a' ,a 

= U6(l, 2)G B (1, l+)<w(l, 1) + V5(l, 2) J2 G °' (1 + a, 1 + a+)g aa , (1, 1 + a) (9) 

<x',a 

where gaa'ihj) is the equal-time pair correlation function which is related to the probability of finding one electron 
with spin a' on site j when another electron with spin a is held on site i. More specifically, 



[ ( _ r r 4(l)c CT (l)4>(2)c CT ,(2))-J(l,2)^(T T 4(l)c CT (l)) K (i)^( 2 ))-J(l,2)^(n CT (l)) 

<?crcr' (,-L, _ . . . . _ (1UJ 



r T 4(l)c CT (l))(T T 4/(2)c^(2) 



(n ff (1)) (n^ (2)) 



r 



where (n a (l)n a ' (2)) is the density-density correlation 
function. In this last formula it is assumed that 
t\ = T2- With this procedure, the four point func- 
tion (T T ct (l)c 5 -(l)c cr (l)4(3)^> appearing in the defini- 
tion of the self-energy Eq.© is factorized a la Hartree- 
Fock everywhere^ except when the point 3 is equal to 1 + , 
in which case there is no approximation involved. The 
Fock contribution from the V term is discussed in Ap- 
pendix A. It gives a very small contribution in the regime 
studied in the present paper. We caution the reader that 
the above-mentioned factorization is not exactly the one 



which is used in the STLS approximation. Further fac- 
torizations and additional details will be discussed in the 
following sections. 

We want to calculate the spin and charge response 
functions. These can be obtained from the first func- 
tional derivative of the Green function with respect to the 
external source field. Taking the functional derivative on 
both sides of the identity G CT (1, 3)G-/(3, 2) = 6 aa ,5(l, 2) 
and using the Dyson equation Eq. 0, we obtain the ex- 
act result 



IW (1,2; 3, 3) = - 



<5G CT (1,2) 
*M3,3) 



(11) 
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In turn, the first functional derivative of the self-energy respect to the Green function can be evaluated from our 
approximate expression for the self-energy Eq. @ 



<SS g (4,5) 
<5G CT »(6,7) 



= [/<W'<5(4, 5)5(4, 6)6(5, 7). 9 ^(4, 4) + V ^ (5(4, 5)5(4 + a, 6)5(5 + a, 7) w (4, 4 + a) 

a 

+ CT( 4,5) G ,(4,4 + )^|i) 4-^(4,5) £ ^(M*) ^**'' - (12, 



r 



The functional derivative of the pair correlation function 
with respect to the Green function is a three-body (six- 
point) correlation function that is not known. For the 
standard Hubbard model, it was shown that the unknown 
functional derivative (third term in the above equation) 
appears only in the charge response function. The au- 
thors in Ref. l(ll l3l|) approximated this functional by a 
constant whose value was obtained by enforcing the Pauli 
principle expressed as a sum-rule on spin and charge cor- 
relation functions. In our case, two other unknown func- 
tional that come from the last term in the above equa- 
tion, appear in both the charge and spin response func- 
tions. We assume, and confirm with the numerical results 



of the following section, that these two unknown func- 
tional do not give important contributions as long as 4V 
is small compared to the bandwidth. Their contribution 
becomes more significant as V increases. By approximat- 
ing the two unknown functions by two different constants, 
it should be possible to obtain them by using two extra 
sum-rules, such as the compressibility and spin suscepti- 
bility sum-rules. We leave this for future work and, at 
this point, we simply drop the <5g CTCT <<' (4, 4 + a)/<5G (T »(6, 7) 
term in the last line. 

The spin and charge part of II can now be obtained by 
combining Eqs. Hll|) and 112fl in the form 



n cc , ss (l, 2; 3, 3) =^(acr')IW (1, 2; 3, 3) = 2 [n ff(T (l, 2; 3, 3) ± Fw(l, 2; 3, 3)] 



2(3,(1, 3)G CT (3, 2) - 2£/G ff (1, 4)G CT (4, 2)^(4, 4) 



5G a (4,4+) SG a (A,4+) 
5^(3,3) <50a(3,3) 



2VG a (l, 4)G CT (4, 2) J2 W(4, 4 + a) 

a' ,a 

2C/G CT (l,4)G ff (4,2)G^(4, 4) 



5G a {l + a~A + a+) £G CT (4 + a,4 + a+)~ 



<%(3,3) 



± 



tya(3,3) 



5.9^(4,4) ± (5^(4,4) 



ff (3,3) 5<M3,3) 



(13) 



where a — ± (the third identity in the above equation is 
valid for a spin-unpolarized system). 

To obtain the response functions of the system, we 
set the external potential to zero. When the minus sign 
(corresponding to the spin response function) is chosen 
in the last term, it drops out by rotational invariance in 
zero source field^. For the plus sign (corresponding to 
the charge response function), we assume that the func- 
tional derivative of the pair correlation respect to the 
density is a constant (after using an extra chain rule in 
the above equation). The final form of the charge and 
spin response functions, or equivalently susceptibilities, 
in Fourier space then have the following forms 



x°(q.^) 



U ss (q) 



(14) 



Xcc(q,^™) 



x°(q,w„ 



1 , x (q-^) 



u c M) 



(15) 



where 



U ss (q) = Ug aa (0) - 4Vg ss {a) J2 cos(q a a), (16) 



C/ cc (q) = U (gaa{0) + n<?U(0)) + 4Vg cc (a) ^ cos(q a a) 

a 

and uj n = 2imT is the Matsubara frequency. One can 
easily see the similarity of the above equations with what 
we had in Refsi^ and also find out how easily these 
equations can be extended to a system with longer range 
of interaction and also to a system with many bands. 
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Finally the index a takes the values a = 1..D, D being 
the dimension of the system, while the equal-time charge 
and spin pair correlation functions are defined by g cc . ss = 
Y,aa>( aa ') n cr n <T'9<7<T' /n? (or simply 5 CC:SS (r) = (^(r) ± 
g a a{v))/2 for a spin-unpolarized system) with |r| = or a 
and x°(q, u) n ) is the response function for non- interacting 
electrons given by 

,o (q ,„„ )= / jp /°(P + f)-/>-?) 

JBZ v tuJ n — e p+q/2 + Cp-q/2 

In the above formula v is the volume of the Brillouin 
zone (BZ), /°(q) = 1/[1 + exp((e g - (J,)/T)\ is the Fermi- 
Dirac distribution function, and e q = —2tJ2 a cos(q a a) is 
the non-interacting particle dispersion relation. The pair 
correlation functions are related to the static structure 
functions by 

5cc(r) = 1 + - [ ^[S cc (q) - 1] exp(jq • r) (19) 
nJ B z v 

g ss (r) = - [ ^ [S ss (q) - 1] exp( iq • r) (20) 
n Jbz v 

where S CCiSS = J2aa' ( aa ' ) V n ^ n ^' &<r<r' l n ( or simply 
S C c,ss(r) = S aa (r) ± SW^r) for a spin-unpolarized sys- 
tem) are respectively the spin and charge static structure 
factors, with the spin-resolved static structure factor de- 
fined by S a(T '(r) = [(7v(0)7v(r)) I y/n a n a >\ - ^n a n a : 
and the Fourier transforms by, 

£e*P- r >^(p), f ^ r 'E 6 ru0 . (21) 
j Jbz v 

Finally the static structure factors are connected to the 
response functions by the fluctuation-dissipation theorem 

S 00 , (q) = _L_ V X**> (q, (22) 



The above equations [Eqs. (|15|I - H22() ] form eight rela- 
tions containing nine unknowns. The extra unknown can 
be fixed using the Pauli principle, namely g CTCr (0) = or 
9cc{0) = ~3ss(0). 

To conclude, note that the RPA approximation on the 
nearest-neighbor interaction V can be simply recovered 
by setting <? cc (l) = 1 and g ss (l) = which means that 
the RPA does not give any extra correction to the spin 
response function of the extended Hubbard model. This 
is a consequence of the fact that the different spin com- 
ponents (spin parallel and anti-parallel) of the off-site 
interaction are identical in the original Hamiltonian. We 
will see that this is in contradiction with QMC calcula- 
tions. Thus, our generalization of TPSC does take into 
account important non-perturbative corrections. Finally, 
at this level of approximation the self-energy Eq. is 
a constant. As in the original TPSC approachi2Lii&, we 
can perform a second step to improve our approximation 
for the self-energy. This is discussed in Appendix 1X1 
B. Wigner distribution function approach 

In this section, we present the approach for obtaining 
the structure factors that parallels that of STLS 8 . It is 
based on the equation of motion for the Wigner distri- 
bution function. We just show the calculation for the 
ordinary Hubbard model (V = 0) in order to shorten the 
length of the equations. This will suffice to demonstrate 
the difference between TPSC and STLS. 

The one- and two-body Wigner distribution function 
(OBWDF and TBWDF) are defined by, 



ha (p, r) = £ e* 1 " (4 +1/2iCT Ci_i/ 2 , CT ) (23) 
l 



I 

/iiw(p,p',T) =£ e *P-i e V "v ( C I +1/2 , CTC i-i/ 2 ,.4 + iV2, CT ' c i-i'/2,.') • (24) 



One should notice that the operators in the above equa- 
tions act on some lattice sites which do not exist in the 
real system. Knowing that the Winger distribution func- 
tions (WDFs) are not a real physical functions we define 
them in this manner for the sake of simplicity in the 
notation. One can also define the WDFs in term of the 
operators which just act on the real lattice sites, but that 
makes the formalism a bit more tedious. 

The density of particles of spin a at position i is related 



to the OBWDF by 

n-Ur) = f — / iff (p,r) (25) 

JBZ v 

with the same definitions of Fourier transforms as above, 
Eq. (HJ. 

We first need to write the equation of motion for the 
operator c\^ a to obtain the equation of motion for the 
OBWDF later on. This equation, after a bit of algebra, 
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can be written as, 

(26) 



where £fci l<T = J2q). c j,<? an d means that the sum 
runs over all nearest-neighbors of site i. Using the above 
equation, one can write the equation of motion for the 
OBWDF, 



<;/i " T(p - r ) = -2i*5^sm(g a )(*i) a /i jCr (p > T) 



Or 



-[/W — / — e l(p - pi)ri (5i-iM/ 2 -5i-i'.-i/ 2 )/iiw(pi,P2,r). (27) 
j Jbz v Jbz v 

To derive the first term on the right-hand side of the above equation, we used the following identity, that we prove in 
one-dimension only (for the sake of simplicity of notation): 



= C i+ ; /2+1)0 .Cj-V2,(7 - C J+Z/2, CT C -! ; -'/2-l,cr + C i + l/ 2 -l, a C i-l/2,c ~ 4+1 / 2 ,o C i~l 1 1+\ ,<J 
= 6 i( C ]+l/2+l/2^-l/2-l/2,a) ~ ^i{c i+ l/ 2 -X/2,<T C i-lft+l/^) 

= Sifal(4 +l/2>a .c i _ l/2ia ), (28) 

where Sihi t<r = h i+ i/ 2a — ^i-1/2,0- and §2%hi. a — hi+i fC r — 7Ji-i jCr where hi is a general function of the operators Cj and 
cj, such as hi = 4+j c i+k in which j and k are arbitrary numbers. 

The TBWDF appears in the equation of motion for the OBWDF Eq. I|27|). which means that we have to make 
an approximation in order to obtain a closed set of equations. Proceeding by analogy with the previous section, wc 
factor the TBWDF as follows, 

/iiw(p,p',r) « / iCT (p,T)./W(p',T). 9CT5 (i',i',T) (29) 
« /<r(p)/*(p')s™(0) + fa(p)fva(p',r) ga9 (Q) + /i CT (p,r)/ ? (p')<7™(0) 

+ Up)Mp') E d9 Z^l(r) ) fii ' AT) (30) 

I 



where we define /i CT (p,T) = / CT (p) + /i CT (p,r) with 
/i<r(p 5 T ) the deviation of the OBWDF from its average 
value due to presence of the external potential. In addi- 
tion we assume that the external potential is weak enough 
that we can keep only the first term in the functional Tay- 
lor expansion of g<r&(i', i', r), which means that /i CT (p, t) 
and n^ a i(r) are small. The first and third terms of the 
above equation do not contribute to the final form of 

I 



the equation of motion for the OBWDF, Eq. The 
functional dependence of g^z (i, i', t) on n- 3 ^i (r) again ap- 
pears in the above equation. We use the well-known local 
approximation dg a „{i, i, t) / dnj >(T > (t) — 6i.jdg a z(0)/dn a i 
where n a is the average number of particles per site with 
spin a. The final non-zero contribution from the above 
approximation for the TBWDF to the equation of the 
motion finally takes the following form, 



x^iw(p ) p^T)=/ CT (p)^ ^ (p^^Ko)+/ CT (p)^(pOE^ £i ^ ni ^-'( r )• ( 31 ) 

The exact form of /o-(p) is not known but in first approximation it is reasonable to replace it by the Fermi-Dirac 
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function /"(p). The final equation for /i CT (p, r) in Fourier space can finally be written as 

K - 4^sin(|)sin(p a )]/ CT (q,p, Wn ) = [/°(p + |) - /°(p - |)]^r*(q, Wn ) 

+ tfW ^^(p-PO-n [e if. ri _ e ^f- ri] ^ CT;(q)P)P ^ n) . (32) 

Y Jbz v Jbz v 

Now we can perform an integral over p to obtain an equation for the density n CT (q, 0J n )- The final form of this equation 
is 

n CT (q,w n ) = xl((i,^n)V^ xt (q,uj n ) + U\l (q, ^n)[n g (q, cj n )g^(0) + ^ ^^^^ ^(q,^)] (33) 



where x°(q, is given by, 

dp/5(p + f)-/S(p-f) 



BZ f ZW n + £p+q/2 ^ ^p_ q /2 



(34) 



One can invert the equation for the change in density 
Eq. (|33|l in order to obtain the density in terms of the 
external potential to extract the susceptibility, 

n ff (q,w„) = Xcc(q,Wn)V c e c a:t (q,w n )+Xs S (q,w„)K 3 e s ;l:t (q,w n ) 

where V c e c %(q,u n ) = [V a ext (q, u n ) ± V£ xt (q,u n )}/2, and 
the coefficients of the external potential are the response 
functions, which are given by following formulae, 



s (q, w„) 



x°(q,^«) 



it 



i x°(q ! Wn) 



(35) 



where X°(<l,^n) = 2 (x°(q,o;„) +X°(q,w n )), C/ ss = 
C/.g CT a(0) and C/ cc = C/[<w(0) + nfl ffcrff (O)/0n]. 

Eqs. I|35|l are the same as Eqs. ()15fl and Eq. (|14|) when 
we set y = 0. The extension of the above equations to 
the case V ^ is straightforward and leads to exactly 
the same result as in the previous section. In the present 
case, the derivative of the pair correlation function with 
respect to the density can be evaluated if one wishes, 
but its contribution is not big enough to reproduce the 
QMC results as we will show in the next section. This 
problem is known in the context of the electron liquid 36 . 
The authors add a extra unknown multiplier constant 
and fix it by the compressibility sum-rule. If they had 
used instead Pauli principle they would have recovered 
the TPSC equations. 



C. Comments on the STLS approximation 

We are now in a position to contrast the results of 
the above section with the STLS approximation. For the 
sake of simplicity, it is preferable to limit ourselves to the 
case V = 0. The factorization of the TBWDF that leads 
to the STLS approximation is given by 

h><Ts(p, P', t) ps / io .(p, T)fvs(p', T)g<rz(i, i', r). (36) 
This should be contrasted with the TPSC factorization 
appearing in Eq. I|29|l where the pair correlation function 



is taken "on-site" . At first glance the STLS factorization 
looks more reasonable because, as far as the TBWDF is 
concerned, the integral of the last formula with respect 
to p and p' leads to the exact result 

{ni a (r)ni'a'(T)) = {ni a (r)) (n Va ,(T)} g a a{h i', r). (37) 

However, one must recall that in the equations of motion, 
the TBWDF appears weighted by the range-dependent 
potential appearing in the Hamiltonian. In particular, 
the form Jw oo (p 5 p'j T ) is valid only for interactions with 
a finite range. With a local interaction, three of the 
creation-annihilation operators are at the same point, as 
can be seen in Eq. (|26|) . The factorization that appears 
correct, as judged by comparisons with QMC, is the one 
that takes the role of the potential into account. In the 
case of the simple Hubbard model, the potential is local 
in time and space so one needs a local factorization to 
model the interaction terms a best as possible. 

The formal STLS approximation can be obtained by 
replacing the above STLS factorization Eq. I|36ll in the 
equation of motion Eq. I|27f) and then by repeating the 
same steps as above. We must also ignore the functional 
derivative of the pair correlation function with respect 
to the density to recover the simplest result. The final 
forms of the response functions are given by 



i(q,u r , 



x°(q,w ri 



lT|[l- G ffS (q,u;„)]x (q, UJ„ 



(38) 



where G cr g-(q, u> n ) is the local field factor for the qSTLS 
approximation. It can be written as 

G a ir(q,u} n ) = -- / — S^k-q) — ^ — (39) 

n Jbz v Xh\% W «J 

where X°(q, k, u n ) is the inhomogeneous free response 
function 



X^(q,k,cV- ' = ~ 



BZ v w n + e p+q/2 ~ e p-q/2 



(40) 



The local field factor in the STLS approximation can be 
obtain by taking the following limit 37 , 
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G^J iS (q)= lim G CT cf(q, w n 



dk„ . E^m^sin^) 

- o CT( n K qj _ . 2 , q - 



This integral can be simplified using k — > k + q so that the final result appears as q independent, 



r,STLS 



(41) 



(42) 




(0,0) 



(q x .q v ) 




FIG. 1: The static structure factors for different methods 
at U — 8, n — 1 and T = 1 as a function of momentum. 
The upper and lower curves are related to spin and charge 
components respectively. 



FIG. 2: The inverse of the spin part of the static response 
function at specific momenta and V = 0, U = 4, n = 1 as a 
function of temperature. Symbols are QMC results extracted 
from I30L 



III. NUMERICAL RESULTS 

We now present our numerical results obtained from 
Eqs. (|15|) - (|22[) . Unless mentioned otherwise, we use U = 
4 and n = 1 in all the figures, in units where h = ks = 
t = 1. We first present V — results to contrast TPSC 
with other approaches and understand the source of the 
differences, then we move to the more general case. 

A. TPSC, STLS and other approaches for V = 0. 

In Fig. 1 we compare the static structure factors for dif- 
ferent methods at U = 8, n = 1 and T = 1. We take the 
TPSC results represented by the solid lines in Fig. 1 (spin 
on top, charge on the bottom) as our reference. Indeed, 
it was shown in great detail before 1,2 , that the TPSC val- 
ues for the spin and charge structure factors agree very 
closely with QMC calculations that are essentially exact 
within small statistical uncertainties. However, TPSC is 



a weak to intermediate coupling method, so it eventually 
fails for U > 8t. Nevertheless, if one is not too close 
to phase transitions, TPSC results for the spin structure 
factor are still excellent at U = 8t while the results for 
the charge structure factor begin to deviate from QMC 
because of the approximation involved in the evaluation 
of the functional derivative. The STLS results, repre- 
sented by the long-dashed line in Fig. 1, deviate substan- 
tially from TPSC. The inaccuracy of the STLS method 
for the Hubbard model comes from the fact that the po- 
tential is local, so one should use the local pair correla- 
tion function to correct the factorization in the equation 
of motion instead of the non-local factoring used in the 
STLS approach, as discussed in Sec. Ill CI In addition, 
since both charge and spin structure factors for STLS 
are larger than for TPSC, it is clear that STLS does not 
satisfy the Pauli principle (n 2 ) = {n a ) (or g CT(T (0) = 0), 
a key requirement for electrons on a lattice at large den- 
sity, where the probability of having two electrons on the 
same site is large. For the charge structure factor, one 
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can compare two more methods with STLS and TPSC. 
The local approximation (LA), represented by dots, con- 
sists in using for the effective interaction in the charge 
channel Eq. U cc = U[g aa {0) + ndg ad {Q) / dn]. We 

call LAO, represented by dot-dashed line in Fig. 1, the ap- 
proximation that neglects completely ndg a a{fy / dn. The 
difference between STLS and LAO is very small. How- 
ever, the difference between TPSC, LA, STLS (or LAO) 
is relatively large, demonstrating the importance of the 
functional derivative in this range of physical parameters. 
In the language of authors involved with the local approx- 
imation, since LA does not provide a satisfactory result, 
the multiplier factor in front of the derivative of the pair 
correlation function respect to the density is necessary. 
When this unknown multiplier is obtained from the Pauli 
sum-rule g&a(0) — instead of from the compressibility 
sum-rule, one recovers TPSC. 

In TPSC, it suffices to know g a „(Q) to obtain the static 
spin structure factor. Could this quantity be determined 
from the compressibility sum-rule instead of from the 
sum- rule relating g a ^{Q) to the integral of the structure 
factor? To answer this question, we show in Fig. 2 the 
TPSC results for the inverse of the static (co = 0) spin 
response function (susceptibility) as a function of tem- 
perature, again for the V — case, compared with QMC 
results of Ref. 1)3(JK Both short and long wave-length lim- 
its show a linear behavior in the intermediate and high 
temperature regimes, exhibiting a Curie law. The devi- 
ations from the Curie law appear at low temperature in 
both QMC and in TPSC. The agreement between TPSC 
and QMC is much better near wave vector (tt, ir), even 
though the deviations are not large even around (0, 0) . 
Hence, the spin susceptibility sum-rule is also very nearly 
satisfied with this method, meaning that the momentum 
independent correction factor, which is given by <? CT a-(0) 
in the effective interaction, corrects the result properly 
over the entire Brillouin zone 38 . However, to use the 
spin susceptibility sum-rule to fix the constant correction 
factor, one needs an independent way to find the spin- 
susceptibility. Normally, the long wave-length behavior 
of the spin (charge) response function is related to the 
second derivative of the free energy with respect to mag- 
netization m = n-f — rii (density n) which can then be 
computed from the free-energy. In TPSC however, the 
free-energy requires further study 3 ^. In addition, given 
the less accurate results exhibited in Fig. 2 near (0, 0), we 
consider it far more preferable to use the original TPSC 
method where <7cra-(0) is determined by an integral over 
all wave vectors (Eqs. <|20ll and (|22|0 so as to satisfy the 
Pauli principle g aiJ {Q) — 0, which also involves a sum 
over all wave vectors. 
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FIG. 3: The staggered static structure factors as a function 
of temperature for U = 4, n = 1, V = 0, 0.5 and 1. The dots 
are the QMC results of Ref. while lines are our results. 
The upper and lower curves (dots) are related to the spin and 
charge components respectively. 
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FIG. 4: The staggered static charge and spin susceptibilities 
as a function of temperature using the same parameters and 
symbols as in Fig. 3. 



B. QMC vs generalization of TPSC for V ± 

To judge the accuracy of TPSC, we plot in Fig. 3 
the staggered static structure factors as a function of 
temperature, using QMC results of Ref. I|30l) as refer- 



ence. The results for the charge structure factor are all 
smaller than the corresponding results for the spin struc- 
ture factor. Our generalization for TPSC is plotted for 
V = 0, 0.5 and 1 while the QMC results, represented by 
symbols, are for V = 0.5 and 1. The figure clearly shows 
good agreement between our results and QMC in the high 
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0.5 1 1.5 2 2.5 3 3.5 4 (0,0) (it,0) (ji.ji) (0,0) 



T (q x ,q y ) 



FIG. 5: The variation of g a a(0) and g C c,ss(a) as a function of 
temperature for U = 4, n = 1 and V — 0, 0.5, 1. 



temperature region or when V = 0.5. The deviation be- 
comes significant when V = 1 and the temperature is low. 
This suggests that the effect of the functional derivative 
that we ignored in writing down Eas. (|15l) and H14[) be- 
comes important when when 4V is of the order of the 
bandwidth (which is 8t here). Indeed the factor of four 
is necessary to account for the number of neighbors. The 
QMC data shows (not on the figure) that the tendency to 
staggered spin order disappears around V ~ 1.25 while in 
our case it persists to a higher value V ~ 2. This means 
that even though g ss (a) decreases with increasing V, the 
combination Vg ss (a) does not increase fast enough. The 
functional derivative must become important to take this 
effect into account. Note that in RPA, the spin structure 
factor is independent of V. 

Finally, Fig. 4 shows the staggered static charge and 
spin susceptibilities as a function of the temperature for 
the same parameters as the previous figure. All features 
are similar to what was mentioned in Fig. 3, except this 
fact that the correction factors do not correct the struc- 
ture functions and response functions in the same man- 
ner. This is due to this assumption that these functions 
are local in time and space which certainly fails in certain 
region of the parameters even in the case V = 0. 



C. The effect of V 

We finally turn to our main point, a more general 
overview of the effect of the nearest-neighbor interaction 
V over a wide range of parameters. In Fig. 5 we show 
the variation of g a d{Q) and g CC;SS (a) as a function of tem- 
perature for V = 0, 0.5 and 1. We first notice that at 



FIG. 6: The charge and spin components of the structure 
factors at T = 1, n = 1, U = 4, V = -0.5, and 0.5. The 
upper and lower curves correspond respectively to the spin 
and charge components. 

V = 0, both <?ctct(0) and g S s(a) have a sharp decrease 
around T ss 0.3. For g a a{0) this means a decrease in the 
probability for finding two particles at the same place. In 
other words, there is an increase in the size of the local 
moment. The fact that g ss (a) = (ffo-o-(a) ~ 9aaifl))/2 is 
negative means that the probability of finding two elec- 
trons at a distance a with opposite spins is larger than 
finding them there with the same spin. The decrease 
of g S s(a) with temperature indicates a reinforcement of 
this tendency. These results reflect the tendency toward 
antifcrromagnetic order (staggered spin order). Long- 
range spin-density wave order occurs only at zero tem- 
perature, as required by the Mermin- Wagner theorem in 
a two dimensional system, but the decrease in g a a(0) 
and g ss (a) reflects the beginning of the renormalized- 
classical regime where the characteristic spin fluctuation 
frequency becomes smaller than temperature and where 
the antifcrromagnetic correlation length begins to in- 
crease exponentially^. By contrast, in the charge chan- 
nel g C c{o) does not show any strong change in the low 
temperature limit so there is no tendency to charge den- 
sity wave order at these values of V within our present 
approximation even at very low temperature. We observe 
that as we increase the V, the staggered spin fluctuations 
are depressed since g a a(0) and g cc ,ss(a) do not decrease 
sharply. 

How spin and charge fluctuations arc influenced by V 
is best illustrated in Fig. 6 where we show the structure 
factors at T = 1, for V = -0.5, and 0.5. Both func- 
tions show a peak around q x = ir and q y = it, a sign of 
the tendency towards staggered ordering. It is obvious 
from the figure that antiferromagnetic fluctuations are 
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FIG. 7: The crossover temperature as a function of filling 
factor for U = 4, and V = -0.4, -0.2, 0, 0.2, 0.4. Positive 
V reduces the strength of antiferromagnetic fluctuations while 
negative V has the opposite trend. 
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At this temperature, the antiferromagnetic cor- 
relation length begins to increase exponentially. 
We define it as the temperature where the ratio 
Xssilx, Qy, 0)/x°{lx, %, 0) at its peak reaches the value 
100. The rounding (foot) of the curves as the crossover 
temperature vanishes comes from this choice of ratio, 
which corresponds to fixing the correlation length at 
which we consider that we have entered the exponential 
regime related to the existence of long-range order at 
zero temperature. When the ratio is taken as larger, 
the rounding occurs over a narrower range of densities. 
In general there are quantum critical points at zero 
temperature where the exponential regime ends2£. In 
Fig. 7, the crossover temperature and the range of 
fillings where antiferromagnetic fluctuations are large 
both increase or decrease together. It is clear from this 
figure that positive V tends to suppress the staggered 
spin fluctuations while negative V leads to the opposite 
trend. As mentioned above, each site has four nearest- 
neighbor sites, and our method is quite accurate as long 
at |4V| < W. 

The apparent breaks in slope on the curves of the 
last figure correspond to a change from commensurate 
to incommensurate fluctuations. We have been using 
the terms antiferromagnetic and staggered rather loosely. 
Fig. 8 corrects this by presenting the peak position at 
the crossover temperature of the spin response function 
(or spin structure factor) as a function of filling factor 
for the same parameters as Fig. 7. The maximum in 
the spin fluctuations changes with temperature but, at 
the crossover temperature, it is at a wave vector that is 
commensurate near half-filling and then becomes incom- 
mensurate as we decrease n. This comes mainly from the 
change in the peak position of the non-interacting sus- 
ceptibility, but when V is finite there is an also effect 
that comes from the presence of the cosine functions in 
the effective interaction appearing in the denominator of 
Eq. (|14l) . The range of fillings where commensurate fluc- 
tuations appear is increased when V < and decreased 
when V > 0, concomitant with the tendency to increase 
or decrease the antiferromagnetic crossover temperature 
at half-filling. 



IV. CONCLUSION AND SUMMARY 



FIG. 8: The wave- vector where the maximum of the spin 
response function appears at the crossover temperature, 
as a function of filling factor for U — 4, and V = 
-0.4, -0.2, 0, 0.2, 0.4. 



suppressed with increasing V while the charge fluctua- 
tions are enhanced. A negative value of V reverses the 
trend. For negative V, pairing fluctuations should also 
become important but they are not considered here. 

In Fig. 7 we show the crossover temperature as a func- 
tion of filling factor for V = -0.4, -0.2, 0, 0.2 and 0.4. 



We generalized the TPSC approach to the extended 
Hubbard model that contains nearest-neighbor repulsion 
V in addition to the usual on-site U term. That non- 
perturbative approach, which is a close relative of the 
STLS approach used for the electron gas, is valid in 
the weak to intermediate-coupling limit. The TPSC ap- 
proach is usually studied in the functional derivative for- 
malism and the STLS approach in the Wigner distribu- 
tion function formalism, so we presented derivations in 
both languages to better illustrate the similarities and 
differences between the two approaches. 

To derive either TPSC or STLS, two main approxima- 
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tions are necessary: First we must factor a four-point cor- 
relation function (two-body density matrix) and correct 
the factorization with the pair correlation function; Sec- 
ond, we must treat the functional derivative of the pair 
correlation function with respect to a fictitious external 
potential. The STLS factorization of the pair correlation 
function does not take into account the range of the inter- 
action whereas TPSC does. In particular, at V = 0, the 
STLS approximation involves the pair correlation func- 
tion for all distances and does not enforce the Pauli prin- 
ciple. On the other hand, TPSC involves the calculation 
of g<ra (0) only and, in addition, it enforces the Pauli prin- 
ciple g aa (0) = 0. The Pauli principle gives an additional 
equation that allows an approximate evaluation of the 
functional derivative entering the charge channel when 
one neglects the momentum and frequency dependence 
of that functional derivative. The local approximation 
(where the functional derivative is replaced by a deriva- 
tive with respect to density) is another approach, but it 
is less accurate than TPSC, as judged from comparisons 
with QMC. Analogous comparisons with the spin and 
charge structure factors obtained from QMC also show 
that TPSC is more accurate than STLS. 

When V ^ 0, one uses the same type of factoriza- 
tion, but extra functional derivatives appear in TPSC. 
These extra derivatives cannot be determined from the 
same kind of simple sum-rules used for the V = Hub- 
bard model. Comparisons of numerical results with QMC 
simulations show that for 4V < W, one can neglect the 
extra functional derivatives. In principle, one could ob- 
tain these derivatives within the local approximation or, 
given an independent way to obtain the free energy, by 



enforcing the compressibility and spin susceptibility sum- 
rules. This is left for future work. 

At V = 0, near half filling, there is a crossover tempera- 
ture to the renormalized-classical regime where spin fluc- 
tuations near (71", w) grow exponentially, diverging only at 
zero temperature in agreement with the Mermin- Wagner 
theorem in two dimensions. The effect of V > is to de- 
crease both the crossover temperature and the range of 
dopings over which this crossover occurs. The crossover 
to increasing commensurate (tt, it) fluctuations occurs 
over a narrower range near n = 1 when V > than when 
V = 0. The opposite conclusions hold true for V < 0. 
All these effects are non-perturbative. They are beyond 
the reach of RPA where V does not influence the spin 
fluctuations. Staggered charge fluctuations are enhanced 
by V > and decreased by V < 0. 
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APPENDIX A: OTHER POSSIBILITY FOR 
FACTORIZATION 

Here we would like to return to Eq. JSJ in order to 
mention the Fock-type factorization of the V term, which 
can be introduced as follows 



Z' (T (l,2) = -Vj2 G Ahl + a)6{l + a,2)g xx {l,l + a) (Al) 

a 

where E^(l,2) is the extra contribution to the self-energy in Eq. @ and g xx (a) is a new pair correlation function 
defined as 

^_ (r T ct(i + a ) c<r (i + a ) c t(i) Cg (i)) 

9xx[a) ~ G CT (l,l + a+)G CT (l + a,l+) ' (A2) 
This pair correlation function is related to the one for parallel spins by 

,x_ gaa{a)nl ... 

gxx[a) - G a {l,l + a+)G a {l +«,!+)• { ' 



The Green function that appeares in the definition of the pair correlation function Eq. (|A2|) can be written as 

G CT (l + a,l+) = J ^e iqa n CT (q) = J ^j[cos(q • a) + i sin(q • a)]n (T (q) = n c (a) + in s (a). (A4) 

where n CT (q) is the momentum distributio n (c L a Cqp-). The extra contribution to the approximate self-energy is 
obtained by substituting Eq. (|A4|) and Eq. <|A3|) in Eq. IjAljl . The final form in Fourier space is 

, 2 n c (a) [cos(q x a) + cos(g^a)] + n s (a)[sm(q x a) + sin(g y a)] 

S ff (q) = -2Vn a g aa {a) — . (A5) 

n c (a) z + n s (ay 

I 



We stress that the n s (a) term is zero if we assume that relation) is a symmetric function. In this case, the fi- 
the momentum distribution function (or the dispersion 
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nal contribution of the above equation renormalizes the 
hopping term as follows 

tl = t+ 2Vnl 9 (a) 
n c [a) 

The effect of this renormalization for positive V is merely 
a depression in the structure functions. The effect is op- 
posit for negative V and it can even lead to an instability 
when t' ~ 0. One also should notice that inclusion of 
n s (a) in the self-energy leads to an asymmetric disper- 

I 



^(4,5) 
5G a „ (6,7) 



By inserting the above equation in Eq. after ignor- 
ing the second term which involves the functional deriva- 

I 



One notices that the three point response function ap- 
pears in the above correction for the two body response 
function. It is hard to determine this type of function 
in practice, but its contribution to the equation for the 
two-body response function can be estimated. Eq. 111|) 

I 



sion relation (or an asymmetric momentum distribution 
function) which in return gives a non-zero value for this 
quantity. The presence of a self-consistent asymmetry 
in the momentum distribution function is known as a 
Pomcranchuk instability, which can be related to pres- 
ence of stripes in the system22iS£. 

The extra term in the self-energy Eq. (JAlfl also leads 
to an extra contribution in the response functions. To 
evaluate this extra correction we need 



(A7) 



tive, one gets the following extra correction to the equa- 
tion for the response functions 



(A8) 

I 

for parallel spins is our basic equation for this task. We 
replace 1 — > 3, + a and 3 — ► 1, then multiply both 

sides of the equation by VG a (l, 3)Go-(3+a, l)fl'a;a;(3, 3+a) 
and finally performing a sum over the internal index 3. 
The final result is given by the following equation 



V5aa» 5 ( 4 + °» 5 W 4 ' 6)5(4 + a, 7)g xx (4, 4 + a) 

a 

sr^ f~, fA A . sxfA . K Jg xx (4,4 + a) 

a 

I 



5G a „ (6,7) ' 



X ' aa M = VG a (l,3)G a (3 + a,l) 



SG a (3,3 + a) 
5^(2,2) 



9xx(3, 3 + a). 



X' aa (h 2) - x£ (2) (l, 2) + Ug aa (0)xZ (2) (h 4)**, (4, 2) 
+ V W (a)xZ (2) (1, 4 + a)x^(4 + a, 2) 

a,a' 

where x™ (2 \^i 2 ) is 

X^ (2) (q, w n ) = Vg xx (a)G a (l, 3)G CT (3, 2)G a (\, 3 + a)G CT (3 + a, 2) (A10) 



while IV stands for the first order in V and (2) is a label Eq. (|A9f) which already contains Xct<t(1; 2). Eq. (|A9|) and 
for its equivalent diagram. One can replace x<j<j' (1, 2) in Eq. (|13fl after equating labels 1 and 2 [in Eq. (|13fl ] form a 



14 



set of coupled equations which in principle can be solved. 
Here we are not going to provide the final form of the 
equations and instead simply explain that the correction 
term is small compared to its counterparts. This can 
be easily understood from Eq. (|A9|) because x<«r^(l, 2) 



J 



appears as a multiplicative factor and also as a separate 
factor in all terms, so having an estimate from this term 
tells us a lot about the importance of the correction term. 
This term can be easily evaluated and it has the following 
form in Fourier space 



xZ (2 \^ n ) = 9xx {a) / — 



dk 



dk' U(k + f, 



/.(k-D/^k' + f) 



/<x(k' 



V lUJ n 



(ek+£ 



r 



e k'-# ; 



V{\k-k'\), 



(All) 



where V(q) = V[cos(q x a)+cos(q y a)]. The above function 
is proportional to one of the first order diagrams that ap- 
pear in the perturbation expansion of the response func- 
tions. We evaluated the above equation numerically and 
compared the result with the other first order diagram, 
which is given by 



(A12) 



We checked that Eq. (| All|> is zero at half-filling within 
our numerical precision and its contribution away from 
half- filling is negligible compared to Eq. (|A12() . 



APPENDIX B: IMPROVED SELF-ENERGY 
(SECOND STEP OF THE TPSC 
APPROXIMATION) 

We can improve our approximation for the self-energy 
to include single-particle scattering off low-energy spin 



J 



and charge fluctuations. These processes give momen- 
tum and frequency dependence to the self-energy2i2i2i^. 
This improved self-energy leads to one-particle spectral 
functions that compare extremely well with QMC in the 
case of the usual Hubbard model 3 . In the extended Hub- 
bard model, the improved self-energy that includes the 
effects of longitudinal fluctuations can be obtained as fol- 
lows. We first write the exact result 



E„(l, 2)G CT (2, 3) = - U (^4(1)09(1)^(1)4(3) 

- (tJAI + o)<v(l + aK (1)4(3) 

a' ,a 

5G a (l,3) 



= -u 



V J2 



(i++,i+) 

<5G CT (1,3) 



(T /(l + a++,l + a+) 
After replacing 3 — > 1 + and using rotational symmetry, we have 



-G ff (l,l + )G« r (l,3) 

-G CT ,(l + a,l + a+)G CT (l,3) 



S CT (1, 2)G CT (2, 1+) = Un a n a g a a{0) + 2Vn 2 g cc {a) 

= U[xaa(l, 1) + n a n a ] + 2V[ X cc(l, 1 + a) + n 2 
I 



(Bl) 



(B2) 



(B3) 
(B4) 



which can be interpreted as a sum-rule relating one par- 
ticle quantities to the left and two-particle quantities to 



the right. Note that the equalities 

g a a{0) = h 1, 

n a n„ 

X ec (l,l + q) 
ffcc(a) = 5 V 1 



(B5) 
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hold. They are the real-space version of the fluctuation- 
dissipation theorem. 

We already approximated the self-energy using a fac- 
torization of Eq. (IB 1|) . Our factorization is exact in the 
limit 3 — * 1 + so we don't expect anything new from that 
equation. However, we can still use an approximate form 
for the response functions in the second identity Eq. (|B2|I . 
Now our main point is that if the response functions sat- 
isfy the sum-rules at the two-particle level Eq. I|B5|I , then 



there is a guarantee that the self-energy satisfies the sum- 
rule in Eq. I|B3[) relating one- and two-particle quantities. 

As an example one can insert Eq. (fill) inside Eq. (|B2|1 . 
which already contains the approximate form of the self- 
energy Eq. 0, to get an approximate form of the self 
energy at a second level of approximation. We give the 
final form of this equation in Fourier space for practical 
use: 



+ C/ cc (q)[J7 + 4^7(q)]xcc(q, w„0}G o (k + q,w„ + u n >). 

I 



(B6) 



where 7(q) = J2 a cos(<7 Q a). As we expect, the above for- 
mula for the self-energy Eq. (|B6fl satisfies the sum-rule 
in Eq. I|B1(I . This, in fact, is a result of using Eqs. (|19|1 . 
(|2U[) . (|21l) and (|22f> . which are another version of the 
fluctuation-dissipation theorem Eq. (|B5|1 . 

We have to mention that we dropped out a term during 
the calculation which is given by 

Khf(1, 2) - -V G a (l, 1 + a)6(l + a, 2). (B7) 



One should include the correction term to the response 
function we already discussed in the last section to take 



care of presence the above term. This means that one 
should include both corrections in the self-energy and the 
response function (or ignore them from both) to have the 
sum-rule Eq. (|B3fl . 



Following the procedure established for the ordinary 
Hubbard model 3-31 , one could also take into account 
transverse spin fluctuations and crossing symmetry to 
write a more general result. 
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